/* Integration points (quadratures)
   Jan Stebel
*/

#ifndef _INTPOINT_H__
#define _INTPOINT_H__

#include <fstream>
#include <string>
#include <iostream>
#include "fem/mesh.h"


/** evaluation of function only, or gradient w.r.t. q or gradient w.r.t. a */
typedef enum
{
  EVAL_F, EVAL_FG, EVAL_FA
} eval_t;


template<unsigned int dim> class Quadrature;
class CStateProblem;





/**
  Definition of Integration (interpolation) points
*/
template<unsigned int dim>
class IntPoint
{
 protected:
  Mesh *imesh;       /**< Associated mesh */
  const Quadrature<dim> *Q; /**< Associated quadrature formula */
  CMat<CVar> J_a,       /**< Jacobian in CVar variables */
             IJ_a,      /**< Inverted jacobian in CVar variables */
	     x_a;       /**< Nodal coordinates in CVar variables */
  CVar IdJ_a,           /**< Determinant of inverted jacobian in CVar variables */
       dJ_a;            /**< Determinant of jacobian in CVar variables */
  CMat<double> SF,      /**< Shape functions */
               GSF,     /**< Gradient of shape functions. */
	       J,       /**< Jacobian */
	       IJ,      /**< Inverted jacobian */
	       x;       /**< Nodal coordinates */
  map<int, map<int, map<int,double> > > HSF; /**< hessian of shape functions */
  double IdJ,           /**< Determinant of inverted jacobian */
         dJ;            /**< Determinant of jacobian */
  CVec dvec;            /**< Auxiliary vector of d.o.f. */
  int *dintvec;         /**< Array of indices of associated d.o.f. */
  int indof;            /**< No. of associated d.o.f. */
//  int dim;              /**< Space dimension of integration point */
  
  /**
    Calculate jacobian.
    @param E Type of evaluation
  */
  virtual void calcJ(eval_t E);
  
  /**
    Evaluate function value and gradient at quadrature point.
    
    @param dofvec Vector of d.o.f.
    @param dof    Array of indices of d.o.f.
    @param val    Returned function value
    @param grad   Returned spatial gradient
    @param E      Flag indicating type of evaluation
  */
  void evaluate(CVec &dofvec, int *dof, CVar &val, CVec &grad, eval_t E);
  
  /** Set gradient of shape functions. */
  virtual void setGSF() = 0;
  
  /** Set shape functions. */
  virtual void setSF()  = 0;
  
  /** True if integration point is linear (i.e. gradient of s.f. is constant) */
  virtual bool isLinear() = 0;
  
  /**
  Set the coordinates of d.o.f.
  @param elem Element number
  @param E    Type of evaluation
   */
  void setX(int elem, eval_t E);

  /** True if the approximating function space is discontinuous */
  virtual bool isDiscontinuous() { return false; }
  
  
 public:
  double ksi[dim]; /**< Local coordinates */
  CVar dx_a;       /**< Integration weight (in CVar type) */
  double dx;       /**< Integration weight */
  
  /** Constructor (not to be used) */
  IntPoint() {};
  
  /**
    Default constructor.
  
    Initializes the integration point and associates it with mesh and quadrature.
    @param mesh Associated mesh.
    @param q    Associated quadrature
  */
  IntPoint(Mesh *mesh, Quadrature<dim> *q) { init(mesh, q); };
  
  /** Destructor */
  virtual ~IntPoint() { delete[] dintvec; };
  
  /** Quadrature getter */
  const Quadrature<dim> *getQ() { return Q; };
  
  /** Getter for dJ */
  double getDJ() { return dJ; };
  
  /**
    Initialization of data structures.
    @param mesh Associated mesh
    @param q    Associated quadrature
    @param ndof Number of d.o.f. per element
    @param dm  Space dimension of integration point
  */
  void         init(Mesh *mesh, const Quadrature<dim> *q, int ndof);
  
  /**
    Virtual method for initialization of particular integration point types.
    @param mesh Associated mesh.
    @param q    Associated quadrature
   */
  virtual void init(Mesh *mesh, const Quadrature<dim> *q) {};
  
  /** Return number of d.o.f. per element. */
  int          getNdof() { return indof; }
  
  /**
    Return number of d.o.f. per element in case that mesh
    consists of different element types.
    @param elem Element number
  */
  virtual int  getNdof(int elem) { return indof; }
  
  /** Return number of mesh nodes. */
  virtual int  getMeshNdof() { return imesh->getNnodes(); }
  
  /**
    Virtual method for calculation of element data.
    @param elem Element number
    @param E    Type of evaluation
  */
  void calcE(int elem, eval_t E);
  
  /**
    Virtual method for calculation of quadrature point data.
    @param qp Quadrature point
    @param E  Type of evaluation
  */
  void calcQP(int qp, eval_t E);
  
  /**
    Virtual method for calculation of quadrature point data at arbitrary point.
    @param l1 Local (barycentric) coordinate
    @param l2 Local (barycentric) coordinate
    @param l3 Local (barycentric) coordinate
  */
  void calcQP(double l1, double l2, double l3, double w, eval_t E);
  
  /**
    Evaluate base function and spatial gradient in quadrature point.
    @param basis  Number of base function
    @param elem   Element number
    @param offset Offset of current variable within global vector of d.o.f.
    @param val    Returned function value
    @param grad   Returned spatial gradient
    @param E      Type of evaluation
  */
  void         evalBasis(int basis, CVar &val, CVec &grad, eval_t E);
  
  void         evalBasisDG(int elem, int parent, int basis, CVar &val, CVec &grad, eval_t E);
  
  /**
    Virtual method for evaluation of function and gradient in quadrature point.
    @param dof    Global vector of d.o.f.
    @param elem   Element number
    @param offset Offset of current variable within global vector of d.o.f.
    @param val    Returned function value
    @param grad   Returned spatial gradient
    @param E      Type of evaluation
  */
  virtual void eval  (CVec &dof, int elem, int offset, CVar &val, CVec &grad, eval_t E);
  
  virtual void evalDG(CVec &dof, int elem, int iParent, int offset, CVar &val, CVec &grad, eval_t E);
  
  /** Return global x-coordinate of actual quadrature point */
  CVar getX();
  /** Return global y-coordinate of actual quadrature point */
  CVar getY();
  /** Return global z-coordinate of actual quadrature point */
  CVar getZ();
  
  /**
    Return global number of d.o.f. given by element and local number.
    @param elem Element number
    @param dof  Local number of d.o.f.
  */
  virtual int mapDof (int elem, int dof); // map local to global index of degree of freedom
  virtual int mapDof(int elem, int dof, int iParent);
  
  friend class CStateProblem;
};










#endif
